#ifndef DORMAND_PRINCE_HPP
#define DORMAND_PRINCE_HPP

#include "EmbRK.hpp"

class DormandPrince : public EmbRK {
public:
  DormandPrince(int dim, double dt, double t0, double t1) : EmbRK(dim, dt, t0, t1) {
    m_s = 7;
    m_ph = 5;
    m_pl = 4;

    if (!initializeSolver())
      std::cerr << "Error in ClassicalRK initialization." << std::endl;
  }

  bool initializeSolver() {
        // 初始化系数
        m_A.resize(m_s, m_s);
        m_c.resize(m_s);
        m_bl.resize(m_s); // 4阶方法系数
        m_bh.resize(m_s); // 5阶方法系数
        
        // 设置c向量 (时间节点)
        m_c << 0.0, 0.2, 0.3, 0.8, 8.0/9.0, 1.0, 1.0;
        
        // 设置A矩阵 (Butcher表)
        m_A.setZero();
        
        // 第1行全为0
        
        // 第2行
        m_A(1, 0) = 1.0/5.0;
        
        // 第3行
        m_A(2, 0) = 3.0/40.0;
        m_A(2, 1) = 9.0/40.0;
        
        // 第4行
        m_A(3, 0) = 44.0/45.0;
        m_A(3, 1) = -56.0/15.0;
        m_A(3, 2) = 32.0/9.0;
        
        // 第5行
        m_A(4, 0) = 19372.0/6561.0;
        m_A(4, 1) = -25360.0/2187.0;
        m_A(4, 2) = 64448.0/6561.0;
        m_A(4, 3) = -212.0/729.0;
        
        // 第6行
        m_A(5, 0) = 9017.0/3168.0;
        m_A(5, 1) = -355.0/33.0;
        m_A(5, 2) = 46732.0/5247.0;
        m_A(5, 3) = 49.0/176.0;
        m_A(5, 4) = -5103.0/18656.0;
        
        // 第7行
        m_A(6, 0) = 35.0/384.0;
        m_A(6, 1) = 0.0;
        m_A(6, 2) = 500.0/1113.0;
        m_A(6, 3) = 125.0/192.0;
        m_A(6, 4) = -2187.0/6784.0;
        m_A(6, 5) = 11.0/84.0;
        
        m_bl(0) = 35.0/384.0;
        m_bl(1) = 0.0;
        m_bl(2) = 500.0/1113.0;
        m_bl(3) = 125.0/192.0;
        m_bl(4) = -2187.0/6784.0;
        m_bl(5) = 11.0/84.0;
        m_bl(6) = 0.0;
        
        m_bh(0) = 5179.0/57600.0;
        m_bh(1) = 0.0;
        m_bh(2) = 7571.0/16695.0;
        m_bh(3) = 393.0/640.0;
        m_bh(4) = -92097.0/339200.0;
        m_bh(5) = 187.0/2100.0;
        m_bh(6) = 1.0/40.0;
        
        return true;
    }
};

#endif